#################################################################################
### Racialized Misinformation, Factual Corrections, and Prejudicial Attitudes ###
### [Main Replication Code for Study 1]                                       ###
### Authors: Eddy S. F. Yeung, Joseph Glasgow                                 ###
### Date: June 26, 2025                                                       ###
#################################################################################

### Set-up ----
## Clean the working environment and set the working directory
rm(list = ls())
setwd("~/Desktop/racialized_misinfo/replication/Study 1") # set your working directory here, which should also contain the survey data ("study1_dataset.csv")

## Import the dataset
df <- read.csv("study1_dataset.csv")

## Load the required packages
library(tidyverse)
library(texreg)
library(estimatr)
library(cowplot)
library(grid)
library(gridExtra)
library(interflex)
library(multcomp)
library(broom)
library(ivreg)
library(modelsummary)
library(effectsize)

### Code demographic variables ----
## Age
df$age <- df$age + 18

## Female (= 1)
df$female <- ifelse(df$sex == 1, 1, 0)

## White (= 1)
df$white <- ifelse(df$race == 1, 1, 0)

## Married (= 1)
df$married <- ifelse(df$marital == 2, 1, 0)

## Education (less than high school = 0; advanced degree = 5)
df <- df %>% mutate(education = case_when(
  educ == 1 ~ 0,
  educ == 2 ~ 1,
  educ == 3 ~ 2,
  educ == 4 ~ 3,
  educ == 5 ~ 4,
  educ == 6 | educ == 7 | educ == 8 ~ 5
))

## Income (less than $10,000 = 0; $150,000 or more = 12)
df$income <- df$income - 1

## Party identification (0 = strong Democrat; 6 = strong Republican)
df <- df %>% mutate(pid = case_when(
  pid_1 == 2 & pid_2d == 1 ~ 0,
  pid_1 == 2 & pid_2d == 2 ~ 1,
  (pid_1 == 3 | pid_1 == 4) & pid_2i == 2 ~ 2,
  (pid_1 == 3 | pid_1 == 4) & pid_2i == 3 ~ 3,
  (pid_1 == 3 | pid_1 == 4) & pid_2i == 1 ~ 4,
  pid_1 == 1 & pid_2r == 2 ~ 5,
  pid_1 == 1 & pid_2r == 1 ~ 6
))

## Democrat (= 1)
df$dem <- ifelse(df$pid == 0 | df$pid == 1 | df$pid == 2, 1, 0)

## Republican (= 1)
df$gop <- ifelse(df$pid == 4 | df$pid == 5 | df$pid == 6, 1, 0)

## Attentiveness (0 = lowest; 2 = highest)
df$screener_2[is.na(df$screener_2)] <- 0
df$screener_2_correct <- ifelse(df$screener_2 == 3, 1, 0)
df$screener_4[is.na(df$screener_4)] <- 0
df$screener_4_correct <- ifelse(df$screener_4 == 2, 1, 0)
df$attentive <- df$screener_2_correct + df$screener_4_correct

## Impute preregistered covariates
df$age[is.na(df$age)] <- mean(df$age, na.rm = TRUE)
df$female[is.na(df$female)] <- mean(df$female, na.rm = TRUE)
df$white[is.na(df$white)] <- mean(df$white, na.rm = TRUE)
df$married[is.na(df$married)] <- mean(df$married, na.rm = TRUE)
df$education[is.na(df$education)] <- mean(df$education, na.rm = TRUE)
df$income[is.na(df$income)] <- mean(df$income, na.rm = TRUE)
df$ideo[is.na(df$ideo)] <- mean(df$ideo, na.rm = TRUE)

### Code outcome variables ----
## Affective polarization (0 = lowest; 1 = highest)
rescale <- function(x) {(x - (-100)) / (100 - (-100))}
df <- df %>% mutate(aff_pol = case_when(
  dem == 1 & group == 0 ~ rescale(affect_C_2 - affect_C_1),
  dem == 1 & group == 1 ~ rescale(affect_T_2 - affect_T_1),
  gop == 1 & group == 0 ~ rescale(affect_C_1 - affect_C_2),
  gop == 1 & group == 1 ~ rescale(affect_T_1 - affect_T_2)
))

## Outparty affect (0 = coldest; 1 = warmest)
df <- df %>% mutate(out_affect = case_when(
  dem == 1 & group == 0 ~ affect_C_1 / 100,
  dem == 1 & group == 1 ~ affect_T_1 / 100,
  gop == 1 & group == 0 ~ affect_C_2 / 100,
  gop == 1 & group == 1 ~ affect_T_2 / 100
))

## Negative stereotype of Blacks in work ethics (0 = weakest; 6 = strongest)
df <- df %>% mutate(stereo = case_when(
  group == 0 ~ stereo_C,
  group == 1 ~ stereo_T
))

## Racial resentment (mean effects index with imputation)
df$rr_1 <- ifelse(df$group == 0, df$rr_C_1, df$rr_T_1)
df$rr_2 <- ifelse(df$group == 0, df$rr_C_2, df$rr_T_2)
df$rr_3 <- ifelse(df$group == 0, df$rr_C_3, df$rr_T_3)
df$rr_4 <- ifelse(df$group == 0, df$rr_C_4, df$rr_T_4)
df$rr_5 <- ifelse(df$group == 0, df$rr_C_5, df$rr_T_5)
rr_matrix <- cbind(df$rr_1, df$rr_2, df$rr_3, df$rr_4, df$rr_5)
cal_MEI <- 
  function(Z, outcome_mat, to_reorient, reorient = F, greedy = T, impute = F) {
    if(impute == T) {
      R <- 1 * is.na(outcome_mat)
      means_for_imputation <- 
        rbind(apply(outcome_mat[Z == 0, ], MAR = 2, FUN = mean, na.rm = T),
              apply(outcome_mat[Z == 1, ], MAR = 2, FUN = mean, na.rm = T))
      to_impute <- R * means_for_imputation[Z + 1, ]
      outcome_mat[is.na(outcome_mat)] <- 0
      outcome_mat <- outcome_mat + to_impute
    }
    c_mean <- apply(X = outcome_mat[Z == 0, ], MARGIN = 2, FUN = mean, na.rm = T)
    c_sd <- apply(X = outcome_mat[Z == 0, ], MARGIN = 2, FUN = sd, na.rm = T)
    z_score <- t(t(sweep(outcome_mat, 2, c_mean)) / c_sd)
    index_numerator <- rowSums(z_score)
    if(greedy == T) {
      n_outcomes <- rowSums(!is.na(z_score))
    }
    else if(greedy == F){
      n_outcomes <- ncol(outcome_mat)
    }
    index <- index_numerator / n_outcomes
    index <- (index - mean(index[Z == 0], na.rm = T)) / sd(index[Z == 0], na.rm = T)
    return(index)
  }
df <- df %>% 
  mutate(resent = cal_MEI(Z = df$group, outcome_mat = rr_matrix, impute = T))

## Attitude toward SNAP (0 = strongly against; 6 = strongly in favor)
df <- df %>% mutate(SNAP_support = case_when(
  SNAP_dec_C == 1 | SNAP_dec_T == 1 ~ 0,
  SNAP_dec_C == 2 | SNAP_dec_T == 2 ~ 1,
  SNAP_keep_C == 2 | SNAP_keep_T == 2 ~ 2,
  SNAP_keep_C == 3 | SNAP_keep_T == 3 ~ 3,
  SNAP_keep_C == 1 | SNAP_keep_T == 1 ~ 4,
  SNAP_inc_C == 2 | SNAP_inc_T == 2 ~ 5,
  SNAP_inc_C == 1 | SNAP_inc_T == 1 ~ 6
))

## Posttreatment belief of the share of SNAP participants (1 = correct)
df <- df %>% mutate(mp_check = case_when(
  group == 0 ~ mp_check_C,
  group == 1 ~ mp_check_T
))

## Posttreatment belief about the percentage of Whites in the US population
df <- df %>% mutate(pop_white = case_when(
  group == 0 ~ population_C_1,
  group == 1 ~ population_T_1
))

## Posttreatment belief about the percentage of Blacks in the US population
df <- df %>% mutate(pop_black = case_when(
  group == 0 ~ population_C_2,
  group == 1 ~ population_T_2
))

## Posttreatment belief about the proportion of SNAP participants among Blacks
df <- df %>% mutate(SNAP_black = case_when(
  group == 0 ~ SNAP_black_C_1,
  group == 1 ~ SNAP_black_T_1
))

## Perceived credibility of the vignette (0 = lowest; 4 = highest)
df <- df %>% mutate(cred = case_when(
  group == 0 ~ cred_C,
  group == 1 ~ cred_T
))

## Overestimation of the share of Black recipients of SNAP (truth = 30.8%)
df$overest <- df$baseline_2 - 30.8

df$group <- as.factor(df$group)
df$treatment <- df$group

### Figure 1 ----
## Average level of misperception in the full sample
mean(df$baseline_2, na.rm = T)
median(df$baseline_2, na.rm = T)
sum(!is.na(df$baseline_2))

## Pretreatment estimate of the proportion of Black SNAP recipients by party ID
group_mean_temp1 <- df %>% 
  filter(gop == 0 | gop == 1) %>% 
  group_by(gop) %>% 
  summarize(xvalue = mean(baseline_2, na.rm = T))
p0.1 <- 
  ggplot(data = subset(df, gop == 0 | gop == 1), 
         aes(x = baseline_2, group = as.factor(gop), fill = as.factor(gop))) +
  geom_density(alpha = .25) +
  geom_vline(data = group_mean_temp1,
             aes(xintercept = xvalue, color = as.factor(gop)), 
             linetype = "dashed", linewidth = 0.5) +
  scale_color_manual(labels = c("Non-GOP", "GOP"), values = c("#00008B", "#EFC000FF")) +
  scale_fill_manual(labels = c("Non-GOP", "GOP"), values = c("#00008B", "#EFC000FF")) +
  geom_vline(xintercept = 30.8, color = "red", linetype = "solid", linewidth = 1) +
  xlab("Estimated % of Black Recipients") +
  ylab("Density") +
  ggtitle("Prior Beliefs by Partisanship") +
  theme_classic() +
  theme(text = element_text(family = "Times", size = 11),
        axis.text = element_text(color = "black", size = 11),
        legend.justification = c(1, 1), legend.position = c(.95, .99),
        legend.box.background = element_rect(color = "black"), 
        legend.key.size = unit(.75, "line"),
        legend.direction = "vertical",
        legend.margin = margin(t = 0.2, r = 0.2, b = 0.2, l = 0.2, unit = "cm"),
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 13, face = "italic"))

## Pretreatment estimate of the proportion of Black SNAP recipients by racial identity
group_mean_temp2 <- df %>% 
  filter(white == 0 | white == 1) %>% 
  group_by(white) %>% 
  summarize(xvalue = mean(baseline_2, na.rm = T))
p0.2 <- 
  ggplot(data = subset(df, white == 0 | white == 1), 
         aes(x = baseline_2, group = as.factor(white), fill = as.factor(white))) +
  geom_density(alpha = .25) +
  geom_vline(data = group_mean_temp2,
             aes(xintercept = xvalue, color = as.factor(white)), 
             linetype = "dashed", linewidth = 0.5) +
  scale_color_manual(labels = c("Nonwhite", "White"), values = c("#00008B", "#EFC000FF")) +
  scale_fill_manual(labels = c("Nonwhite", "White"), values = c("#00008B", "#EFC000FF")) +
  geom_vline(xintercept = 30.8, color = "red", linetype = "solid", linewidth = 1) +
  xlab("Estimated % of Black Recipients") +
  ylab("Density") +
  ggtitle("Prior Beliefs by Racial Identity") +
  theme_classic() +
  theme(text = element_text(family = "Times", size = 11),
        axis.text = element_text(color = "black", size = 11),
        legend.justification = c(1, 1), legend.position = c(.95, .99),
        legend.box.background = element_rect(color = "black"), 
        legend.key.size = unit(.75, "line"),
        legend.direction = "vertical",
        legend.margin = margin(t = 0.2, r = 0.2, b = 0.2, l = 0.2, unit = "cm"),
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 13, face = "italic"))
prior <- plot_grid(p0.1, p0.2, labels = "AUTO", nrow = 1, label_fontfamily = "Times")
ggsave(file = "Figure 1.pdf", prior, width = 8, height = 3)
# ggsave(file = "Figure 1.tif", prior, dpi = 300, width = 8, height = 3)

### Figure 2 and Figure S10 ----
## Work ethic stereotype
# Estimate the ATE
stereo_ATE <- lm_robust(stereo ~ treatment + age + female + white + married +
                          education + income + ideo + pid,
                        data = df)

# Estimate the effect size of the ATE
effectsize::cohens_d(stereo ~ treatment, data = df)

# Compute the general distributions
stereo_by_group <- df %>% 
  filter(stereo >= 0 & stereo <= 6) %>% 
  group_by(treatment, stereo) %>% 
  summarize(count = n()) %>% 
  mutate(freq = formattable::percent(count / sum(count)))

# Visualize the results
group_mean_temp1 <- df %>% 
  group_by(treatment) %>% 
  summarize(xvalue = mean(stereo, na.rm = T))
p1 <- 
  ggplot(data = stereo_by_group, 
         aes(x = stereo, y = freq * 100, group = treatment, fill = treatment)) +
  geom_bar(stat = "identity", position = position_dodge(.9), 
           color = "black", alpha = 0.25) +
  geom_vline(data = group_mean_temp1,
             aes(xintercept = xvalue, color = treatment), 
             linetype = "dashed", linewidth = 0.6) +
  scale_color_manual(labels = c("Control", "Treatment"), values = c("#00008B", "#EFC000FF")) +
  scale_fill_manual(labels = c("Control", "Treatment"), values = c("#00008B", "#EFC000FF")) +
  xlab("Work Ethic Stereotypes against Blacks\n(0 = weakest; 7 = strongest)") +
  ylab("Percentage of Respondents (%)") +
  ggtitle("Work Ethic Stereotypes in Full Sample") +
  theme_classic() +
  theme(text = element_text(family = "Times", size = 13),
        axis.text = element_text(color = "black", size = 13),
        legend.justification = c(1, 1), legend.position = c(.95, .99),
        legend.box.background = element_rect(color = "black"), 
        legend.key.size = unit(1, "line"),
        legend.direction = "vertical",
        legend.margin = margin(t = 0.2, r = 0.2, b = 0.2, l = 0.2, unit = "cm"),
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 13, face = "italic"))
p1 <- p1 + 
  annotate(geom = "text", x = 2.65, y = 33, 
              label = paste("beta == ", round(stereo_ATE$coefficients[2], 2)),
              size = 4, family = "Times", parse = T) +
  annotate(geom = "text", x = 2.65, y = 30.7, 
           label = paste("(SE == ", round(stereo_ATE$std.error[2], 2), ")"),
           size = 4, family = "Times", parse = T)

# HTE by seven-point party ID
out <- interflex(Y = "stereo", D = "treatment", X = "pid",
                 Z = c("age", "female", "white", "married", "education", "income", "ideo"), 
                 data = df, na.rm = TRUE,
                 estimator = "binning", vcov.type = "robust")
p1_HTE <- 
  plot(out, theme.bw = TRUE, cex.axis = 0.8, cex.lab = 0.8, bin.labs = TRUE,
       xlab = "Party Identification\n(0 = strong Democrat; 6 = strong Republican)", 
       ylab = "Marginal Effect on Work Ethic Stereotypes",
       ylim = c(-.75, .75)) +
  ggtitle("HTE on Work Ethic Stereotypes by Party Identification") + 
  theme(plot.title = element_text(hjust = 0.85, size = 13, face = "italic", family = "Times"))

## Racial resentment
# Estimate the ATE
resent_ATE <- lm_robust(resent ~ treatment + age + female + white + married +
                          education + income + ideo + pid,
                        data = df)

# Estimate the effect size of the ATE
effectsize::cohens_d(resent ~ treatment, data = df)

# Visualize the results
group_mean_temp2 <- df %>% 
  group_by(treatment) %>% 
  summarize(xvalue = mean(resent, na.rm = T))
p2 <- 
  ggplot(data = df, aes(x = resent, group = treatment, fill = treatment)) +
  geom_density(alpha = .25) +
  geom_vline(data = group_mean_temp2,
             aes(xintercept = xvalue, color = treatment), 
             linetype = "dashed", linewidth = 0.6) +
  scale_color_manual(labels = c("Control", "Treatment"), values = c("#00008B", "#EFC000FF")) +
  scale_fill_manual(labels = c("Control", "Treatment"), values = c("#00008B", "#EFC000FF")) +
  xlab("Racial Resentment\n(Mean Effects Index)") +
  ylab("Density") +
  ggtitle("Racial Resentment in Full Sample") +
  theme_classic() +
  theme(text = element_text(family = "Times", size = 13),
        axis.text = element_text(color = "black", size = 13),
        legend.justification = c(1, 1), legend.position = c(.95, .99),
        legend.box.background = element_rect(color = "black"), 
        legend.key.size = unit(1, "line"),
        legend.direction = "vertical",
        legend.margin = margin(t = 0.2, r = 0.2, b = 0.2, l = 0.2, unit = "cm"),
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 13, face = "italic"))
p2 <- p2 + 
  annotate(geom = "text", x = -0.60, y = 0.40, 
           label = paste("beta == ", round(resent_ATE$coefficients[2], 2)),
           size = 4, family = "Times", parse = T) +
  annotate(geom = "text", x = -0.60, y = 0.372, 
           label = paste("(SE == ", round(resent_ATE$std.error[2], 2), ")"),
           size = 4, family = "Times", parse = T)
main_ATE <- plot_grid(p1, p2, labels = "AUTO", nrow = 1, label_fontfamily = "Times")
ggsave(file = "Figure 2.pdf", main_ATE, width = 9, height = 4.5)
# ggsave(file = "Figure 2.tif", main_ATE, dpi = 300, width = 9, height = 4.5)

# HTE by seven-point party ID
out <- interflex(Y = "resent", D = "treatment", X = "pid",
                 Z = c("age", "female", "white", "married", "education", "income", "ideo"), 
                 data = df, na.rm = TRUE,
                 estimator = "binning", vcov.type = "robust")
p2_HTE <- 
  plot(out, theme.bw = TRUE, cex.axis = 0.8, cex.lab = 0.8, bin.labs = TRUE,
       xlab = "Party Identification\n(0 = strong Democrat; 6 = strong Republican)", 
       ylab = "Marginal Effect on Racial Resentment",
       ylim = c(-.75, .75)) +
  ggtitle("HTE on Racial Resentment by Party Identification") + 
  theme(plot.title = element_text(hjust = 0.85, size = 13, face = "italic", family = "Times"))
HTE_pid <- plot_grid(p1_HTE, p2_HTE, labels = "AUTO", nrow = 1, label_fontfamily = "Times")
ggsave(file = "Figure S10.pdf", HTE_pid, width = 9, height = 4.5)

### Figure S19 ----
## Support for SNAP
# Estimate the ATE
ATE_SNAP_support <- lm_robust(SNAP_support ~ treatment + age + female + white + 
                                married + education + income + ideo + pid,
                              data = df)

# Estimate the effect size of the ATE
effectsize::cohens_d(SNAP_support ~ treatment, data = df)

# Visualize the ATE
summary_SNAP_support <- df %>% 
  group_by(treatment) %>% 
  do(tidy(lm_robust(SNAP_support ~ 1, data = .))) %>% 
  mutate(SNAP_support = estimate)
p <- 
  ggplot(summary_SNAP_support, aes(x = treatment, y = SNAP_support)) +
  geom_point(data = df, size = 1, alpha = .05, color = "grey",
             position = position_jitter(width = .1, height = .2, seed = 1)) +
  geom_point(size = 1.5) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), linewidth = .5, width = 0) +
  theme_classic() +
  scale_x_discrete(labels = c("0" = "Control", "1" = "Treatment")) +
  xlab("") +
  ylab("Support for SNAP\n(0 = lowest; 6 = highest)") +
  coord_cartesian(ylim = c(-.25, 6.25)) +
  theme(text = element_text(family = "Times", size = 12),
        axis.text = element_text(family = "Times", size = 12))
p <- p + 
  annotate(geom = "text", x = 1.5, y = 5.15, 
           label = paste("beta == ", round(ATE_SNAP_support$coefficients[2], 2)),
           size = 4, family = "Times", parse = T) +
  annotate(geom = "text", x = 1.5, y = 4.85, 
           label = paste("(SE == ", round(ATE_SNAP_support$std.error[2], 2), ")"),
           size = 4, family = "Times", parse = T)
ggsave(file = "Figure S19.pdf", width = 5, height = 5)

### Table S5 ----
## Identify respondents who believed that Black SNAP recipients outnumbered White SNAP recipients
df$misper <- ifelse(df$baseline_2 > df$baseline_1, 1, 0)

## Analyze the correlations between misperceptions and racial stereotype / racial resentment
mod1.1 <- lm_robust(stereo ~ treatment, 
                    data = df)
mod1.2 <- lm_robust(stereo ~ treatment + age + female + white + 
                      married + education + income + ideo + pid, 
                    data = df)
mod1.3 <- lm_robust(stereo ~ treatment + age + female + white + 
                      married + education + income + ideo + pid + misper, 
                    data = df)
mod1.4 <- lm_robust(stereo ~ treatment + age + female + white + 
                      married + education + income + ideo + pid + misper, 
                    data = df, subset = attentive == 2)
mod2.1 <- lm_robust(resent ~ treatment, 
                    data = df)
mod2.2 <- lm_robust(resent ~ treatment + age + female + white + 
                      married + education + income + ideo + pid, 
                    data = df)
mod2.3 <- lm_robust(resent ~ treatment + age + female + white + 
                      married + education + income + ideo + pid + misper, 
                    data = df)
mod2.4 <- lm_robust(resent ~ treatment + age + female + white + 
                      married + education + income + ideo + pid + misper, 
                    data = df, subset = attentive == 2)

## Tabularize the results
texreg(list(mod1.1, mod1.2, mod1.3, mod1.4, mod2.1, mod2.2, mod2.3, mod2.4),
       custom.coef.names = c("Intercept", "Treatment (= 1)", "Age", 
                             "Female (= 1)", "White (= 1)", "Married (= 1)", 
                             "Education (5 = highest)", "Income (12 = highest)", 
                             "Self-Reported Ideology", "Party Identification",
                             "Misperception (= 1)"),
       custom.header = list("Dependent Variable" = 1:8),
       custom.model.names = c("Stereotype", "Stereotype", "Stereotype", "Stereotype",
                              "Resentment", "Resentment", "Resentment", "Resentment"),
       stars = numeric(0),
       include.adjrs = TRUE, include.ci = FALSE, include.rmse = FALSE, 
       caption = "Misperceptions about SNAP Recipiency, Factual Corrections, and Racial Prejudice")

### Figure S11 ----
## Standardize the outcome variables
df <- df %>% mutate(
  stereo_z = scale(stereo),
  resent_z = scale(resent),
  SNAP_support_z = scale(SNAP_support),
  aff_pol_z = scale(aff_pol),
  out_affect_z = scale(out_affect)
)

## Create an empty data frame to store the estimates
df_race <- as.data.frame(matrix(NA, nrow = 10, ncol = 5))
df_race <- df_race %>% 
  rename(race = V1, dv = V2, cate = V3, lower_ci = V4, upper_ci = V5)
df_race$race <- c("Nonwhite", "White", "Nonwhite", "White", "Nonwhite", "White", 
                  "Nonwhite", "White", "Nonwhite", "White")
df_race$dv <- c(1, 1, 2, 2, 3, 3, 4, 4, 5, 5)

## DV = work ethic stereotypes
temp1 <- lm_robust(stereo_z ~ treatment * white + age + female + white + married +
                     education + income + ideo + pid,
                   data = df)
CATE_temp1_nonwhite <- glht(temp1, linfct = c("treatment1 = 0"))
CATE_temp1_white <- glht(temp1, linfct = c("treatment1 + treatment1:white = 0"))
df_race[1, 3:5] <- confint(CATE_temp1_nonwhite)$confint
df_race[2, 3:5] <- confint(CATE_temp1_white)$confint

## DV = racial resentment
temp2 <- lm_robust(resent_z ~ treatment * white + age + female + white + married +
                     education + income + ideo + pid,
                   data = df)
CATE_temp2_nonwhite <- glht(temp2, linfct = c("treatment1 = 0"))
CATE_temp2_white <- glht(temp2, linfct = c("treatment1 + treatment1:white = 0"))
df_race[3, 3:5] <- confint(CATE_temp2_nonwhite)$confint
df_race[4, 3:5] <- confint(CATE_temp2_white)$confint

## DV = SNAP support
temp3 <- lm_robust(SNAP_support_z ~ treatment * white + age + female + white + married +
                     education + income + ideo + pid,
                   data = df)
CATE_temp3_nonwhite <- glht(temp3, linfct = c("treatment1 = 0"))
CATE_temp3_white <- glht(temp3, linfct = c("treatment1 + treatment1:white = 0"))
df_race[5, 3:5] <- confint(CATE_temp3_nonwhite)$confint
df_race[6, 3:5] <- confint(CATE_temp3_white)$confint

## DV = affective polarization
temp4 <- lm_robust(aff_pol_z ~ treatment * white + age + female + white + married +
                     education + income + ideo + pid,
                   data = df)
CATE_temp4_nonwhite <- glht(temp4, linfct = c("treatment1 = 0"))
CATE_temp4_white <- glht(temp4, linfct = c("treatment1 + treatment1:white = 0"))
df_race[7, 3:5] <- confint(CATE_temp4_nonwhite)$confint
df_race[8, 3:5] <- confint(CATE_temp4_white)$confint

## DV = outparty affect
temp5 <- lm_robust(out_affect_z ~ treatment * white + age + female + white + married +
                     education + income + ideo + pid,
                   data = df)
CATE_temp5_nonwhite <- glht(temp5, linfct = c("treatment1 = 0"))
CATE_temp5_white <- glht(temp5, linfct = c("treatment1 + treatment1:white = 0"))
df_race[9, 3:5] <- confint(CATE_temp5_nonwhite)$confint
df_race[10, 3:5] <- confint(CATE_temp5_white)$confint

## Visualize the CATEs by racial identity and by dependent variable (all results)
df_race$dv <- factor(df_race$dv,
                     levels = c(1:5),
                     labels = c("Work Ethic Stereotypes", "Racial Resentment", "SNAP Support",
                                "Affective Polarization", "Outparty Affect"))
p <- 
  ggplot(data = df_race, aes(x = dv, y = cate, color = race, shape = race)) +
  geom_point(position = position_dodge(.25), size = 3) +
  scale_color_manual(values = c("grey25", "grey75")) +
  scale_shape_manual(values = c(16, 17)) +
  geom_errorbar(width = 0, aes(ymin = lower_ci, ymax = upper_ci), 
                position = position_dodge(.25)) +
  xlab("") + 
  ylab("Conditional Average Treatment Effect") +
  theme_classic() +
  theme(text = element_text(family = "Times", size = 13),
        axis.text = element_text(color = "black", size = 13),
        legend.justification = c(1, 1), legend.position = c(.99, .99),
        legend.box.background = element_rect(color = "black"), 
        legend.key.size = unit(1, "line"),
        legend.direction = "vertical",
        legend.margin = margin(t = 0.2, r = 0.2, b = 0.2, l = 0.2, unit = "cm"),
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 13, face = "italic")) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  coord_cartesian(ylim = c(-.5, .5))
ggsave(file = "Figure S11.pdf", width = 10, height = 4)

### Figure S4 ----
## % of respondents disagreeing that "Black recipients > White recipients" in full sample
df %>% group_by(treatment) %>% summarize(n = n(), xvalue = mean(mp_check, na.rm = T))
lm_robust(mp_check ~ treatment + age + female + white + married +
            education + income + ideo + pid,
          data = df)

## % of respondents disagreeing that "Black recipients > White recipients" by party ID
df_manip <- as.data.frame(matrix(NA, nrow = 4, ncol = 5))
df_manip <- df_manip %>% 
  rename(gop = V1, group = V2, percent = V3, lower_ci = V4, upper_ci = V5)
df_manip$gop <- c("Republican", "Republican", "Non-Republican", "Non-Republican")
df_manip$group <- c("Control", "Treatment", "Control", "Treatment")
temp1 <- prop.test(x = sum(df$mp_check[df$treatment == 0 & df$gop == 1], na.rm = T),
                   n = sum(!is.na(df$mp_check[df$treatment == 0 & df$gop == 1])),
                   correct = F)
temp2 <- prop.test(x = sum(df$mp_check[df$treatment == 1 & df$gop == 1], na.rm = T),
                   n = sum(!is.na(df$mp_check[df$treatment == 1 & df$gop == 1])),
                   correct = F)
temp3 <- prop.test(x = sum(df$mp_check[df$treatment == 0 & df$gop == 0], na.rm = T),
                   n = sum(!is.na(df$mp_check[df$treatment == 0 & df$gop == 0])),
                   correct = F)
temp4 <- prop.test(x = sum(df$mp_check[df$treatment == 1 & df$gop == 0], na.rm = T),
                   n = sum(!is.na(df$mp_check[df$treatment == 1 & df$gop == 0])),
                   correct = F)
df_manip[1, 3] <- temp1$estimate
df_manip[2, 3] <- temp2$estimate
df_manip[3, 3] <- temp3$estimate
df_manip[4, 3] <- temp4$estimate
df_manip[1, 4:5] <- temp1$conf.int
df_manip[2, 4:5] <- temp2$conf.int
df_manip[3, 4:5] <- temp3$conf.int
df_manip[4, 4:5] <- temp4$conf.int
df_manip$gop <- factor(df_manip$gop, levels = c("Republican", "Non-Republican"))
p_manip_pid <- 
  ggplot(data = df_manip, aes(x = gop, y = percent * 100, fill = group)) +
  geom_bar(stat = "identity", position = position_dodge(), color = "black", alpha = .5) +
  scale_x_discrete(breaks = c("Republican", "Non-Republican")) +
  scale_fill_manual(values = c("grey15", "grey85")) +
  geom_errorbar(width = .2, aes(ymin = lower_ci * 100, ymax = upper_ci * 100), 
                position = position_dodge(.9)) +
  xlab("") + 
  ylab("Percent of Respondents Disagreeing that\nBlack SNAP Recipients > White SNAP Recipients") +
  ggtitle("Direct Manipulation Check by Partisanship") +
  theme_classic() +
  theme(text = element_text(family = "Times", size = 13),
        axis.text = element_text(color = "black", size = 13),
        legend.justification = c(1, 1), legend.position = c(.34, .99),
        legend.box.background = element_rect(color = "black"), 
        legend.key.size = unit(1, "line"),
        legend.direction = "vertical",
        legend.margin = margin(t = 0.2, r = 0.2, b = 0.2, l = 0.2, unit = "cm"),
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 13, face = "italic")) +
  coord_cartesian(ylim = c(0, 100))

## % of respondents disagreeing that "Black recipients > White recipients" by racial identity
df_manip <- as.data.frame(matrix(NA, nrow = 4, ncol = 5))
df_manip <- df_manip %>% 
  rename(gop = V1, group = V2, percent = V3, lower_ci = V4, upper_ci = V5)
df_manip$gop <- c("White", "White", "Nonwhite", "Nonwhite")
df_manip$group <- c("Control", "Treatment", "Control", "Treatment")
temp1 <- prop.test(x = sum(df$mp_check[df$treatment == 0 & df$white == 1], na.rm = T),
                   n = sum(!is.na(df$mp_check[df$treatment == 0 & df$white == 1])),
                   correct = F)
temp2 <- prop.test(x = sum(df$mp_check[df$treatment == 1 & df$white == 1], na.rm = T),
                   n = sum(!is.na(df$mp_check[df$treatment == 1 & df$white == 1])),
                   correct = F)
temp3 <- prop.test(x = sum(df$mp_check[df$treatment == 0 & df$white == 0], na.rm = T),
                   n = sum(!is.na(df$mp_check[df$treatment == 0 & df$white == 0])),
                   correct = F)
temp4 <- prop.test(x = sum(df$mp_check[df$treatment == 1 & df$white == 0], na.rm = T),
                   n = sum(!is.na(df$mp_check[df$treatment == 1 & df$white == 0])),
                   correct = F)
df_manip[1, 3] <- temp1$estimate
df_manip[2, 3] <- temp2$estimate
df_manip[3, 3] <- temp3$estimate
df_manip[4, 3] <- temp4$estimate
df_manip[1, 4:5] <- temp1$conf.int
df_manip[2, 4:5] <- temp2$conf.int
df_manip[3, 4:5] <- temp3$conf.int
df_manip[4, 4:5] <- temp4$conf.int
df_manip$gop <- factor(df_manip$gop, levels = c("White", "Nonwhite"))
p_manip_race <- 
  ggplot(data = df_manip, aes(x = gop, y = percent * 100, fill = group)) +
  geom_bar(stat = "identity", position = position_dodge(), color = "black", alpha = .5) +
  scale_x_discrete(breaks = c("White", "Nonwhite")) +
  scale_fill_manual(values = c("grey15", "grey85")) +
  geom_errorbar(width = .2, aes(ymin = lower_ci * 100, ymax = upper_ci * 100), 
                position = position_dodge(.9)) +
  xlab("") + 
  ylab("Percent of Respondents Disagreeing that\nBlack SNAP Recipients > White SNAP Recipients") +
  ggtitle("Direct Manipulation Check by Racial Identity") +
  theme_classic() +
  theme(text = element_text(family = "Times", size = 13),
        axis.text = element_text(color = "black", size = 13),
        legend.justification = c(1, 1), legend.position = c(.34, .99),
        legend.box.background = element_rect(color = "black"), 
        legend.key.size = unit(1, "line"),
        legend.direction = "vertical",
        legend.margin = margin(t = 0.2, r = 0.2, b = 0.2, l = 0.2, unit = "cm"),
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 13, face = "italic")) +
  coord_cartesian(ylim = c(0, 100))
manip_check <- plot_grid(p_manip_pid, p_manip_race, labels = "AUTO", nrow = 1, label_fontfamily = "Times")
ggsave(file = "Figure S4.pdf", manip_check, width = 9, height = 4.5)

### Figure S5 ----
## Posttreatment estimate of percent of SNAP recipients among Blacks in full sample
df %>% group_by(treatment) %>% summarize(n = n(), xvalue = mean(SNAP_black, na.rm = T))
lm_robust(SNAP_black ~ treatment + age + female + white + married +
            education + income + ideo + pid,
          data = df)

## Posttreatment estimate of percent of SNAP recipients among Blacks by party ID
summary_SNAP_black <- df %>% 
  group_by(treatment, gop) %>% 
  do(tidy(lm_robust(SNAP_black ~ 1, data = .))) %>% 
  mutate(SNAP_black = estimate) %>% 
  na.omit()
summary_SNAP_black$gop <- factor(summary_SNAP_black$gop, levels = c(1, 0), 
                                 labels = c("Republican", "Non-Republican"))
summary_SNAP_black$treatment <- factor(summary_SNAP_black$treatment, levels = c(0, 1), 
                                       labels = c("Control", "Treatment"))
p_manip_pid_2 <- 
  ggplot(data = summary_SNAP_black, aes(x = gop, y = SNAP_black, fill = treatment)) +
  geom_bar(stat = "identity", position = position_dodge(), color = "black", alpha = .5) +
  scale_x_discrete(breaks = c("Republican", "Non-Republican")) +
  scale_fill_manual(values = c("grey15", "grey85")) +
  geom_errorbar(width = .2, aes(ymin = conf.low, ymax = conf.high), 
                position = position_dodge(.9)) +
  xlab("") + 
  ylab("Estimated % of SNAP recipients\namong the Black Population") +
  ggtitle("Posterior Beliefs by Partisanship") +
  theme_classic() +
  theme(text = element_text(family = "Times", size = 13),
        axis.text = element_text(color = "black", size = 13),
        legend.justification = c(1, 1), legend.position = c(.34, .99),
        legend.box.background = element_rect(color = "black"), 
        legend.key.size = unit(1, "line"),
        legend.direction = "vertical",
        legend.margin = margin(t = 0.2, r = 0.2, b = 0.2, l = 0.2, unit = "cm"),
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 13, face = "italic")) +
  coord_cartesian(ylim = c(0, 100))

## Posttreatment estimate of percent of SNAP recipients among Blacks by racial identity
summary_SNAP_black <- df %>% 
  subset(white == 0 | white == 1) %>% 
  group_by(treatment, white) %>% 
  do(tidy(lm_robust(SNAP_black ~ 1, data = .))) %>% 
  mutate(SNAP_black = estimate) %>% 
  na.omit()
summary_SNAP_black$white <- factor(summary_SNAP_black$white, levels = c(1, 0), 
                                   labels = c("White", "Nonwhite"))
summary_SNAP_black$treatment <- factor(summary_SNAP_black$treatment, levels = c(0, 1), 
                                       labels = c("Control", "Treatment"))
p_manip_race_2 <- 
  ggplot(data = summary_SNAP_black, aes(x = white, y = SNAP_black, fill = treatment)) +
  geom_bar(stat = "identity", position = position_dodge(), color = "black", alpha = .5) +
  scale_x_discrete(breaks = c("White", "Nonwhite")) +
  scale_fill_manual(values = c("grey15", "grey85")) +
  geom_errorbar(width = .2, aes(ymin = conf.low, ymax = conf.high), 
                position = position_dodge(.9)) +
  xlab("") + 
  ylab("Estimated % of SNAP recipients\namong the Black Population") +
  ggtitle("Posterior Beliefs by Racial Identity") +
  theme_classic() +
  theme(text = element_text(family = "Times", size = 13),
        axis.text = element_text(color = "black", size = 13),
        legend.justification = c(1, 1), legend.position = c(.34, .99),
        legend.box.background = element_rect(color = "black"), 
        legend.key.size = unit(1, "line"),
        legend.direction = "vertical",
        legend.margin = margin(t = 0.2, r = 0.2, b = 0.2, l = 0.2, unit = "cm"),
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 13, face = "italic")) +
  coord_cartesian(ylim = c(0, 100))
manip_check <- plot_grid(p_manip_pid_2, p_manip_race_2, labels = "AUTO", nrow = 1, label_fontfamily = "Times")
ggsave(file = "Figure S5.pdf", manip_check, width = 9, height = 4.5)

### Figure S21 ----
## Perceived credibility in full sample
df %>% group_by(treatment) %>% summarize(n = n(), xvalue = mean(cred, na.rm = T))
lm_robust(cred ~ treatment + age + female + white + married +
            education + income + ideo + pid,
          data = df)

## Perceived credibility of the vignette by partisanship
summary_credibility <- df %>% 
  group_by(treatment, gop) %>% 
  do(tidy(lm_robust(cred ~ 1, data = .))) %>% 
  mutate(cred = estimate) %>% 
  na.omit()
summary_credibility$gop <- factor(summary_credibility$gop, levels = c(1, 0), 
                                  labels = c("Republican", "Non-Republican"))
summary_credibility$treatment <- factor(summary_credibility$treatment, levels = c(0, 1), 
                                        labels = c("Control", "Treatment"))
p_cred_pid <- 
  ggplot(data = summary_credibility, aes(x = gop, y = cred, fill = treatment)) +
  geom_bar(stat = "identity", position = position_dodge(), color = "black", alpha = .5) +
  scale_x_discrete(breaks = c("Republican", "Non-Republican")) +
  scale_fill_manual(values = c("grey15", "grey85")) +
  geom_errorbar(width = .2, aes(ymin = conf.low, ymax = conf.high), 
                position = position_dodge(.9)) +
  xlab("") + 
  ylab("Perceived Credibility of Information\n(0 = very untrustworthy; 4 = very trustworthy)") +
  ggtitle("Perceived Credibility by Partisanship") +
  theme_classic() +
  theme(text = element_text(family = "Times", size = 13),
        axis.text = element_text(color = "black", size = 13),
        legend.justification = c(1, 1), legend.position = c(.34, .99),
        legend.box.background = element_rect(color = "black"), 
        legend.key.size = unit(1, "line"),
        legend.direction = "vertical",
        legend.margin = margin(t = 0.2, r = 0.2, b = 0.2, l = 0.2, unit = "cm"),
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 13, face = "italic")) +
  coord_cartesian(ylim = c(0, 4))

## Perceived credibility of the vignette by racial identity
summary_credibility <- df %>% 
  subset(white == 0 | white == 1) %>% 
  group_by(treatment, white) %>% 
  do(tidy(lm_robust(cred ~ 1, data = .))) %>% 
  mutate(cred = estimate) %>% 
  na.omit()
summary_credibility$white <- factor(summary_credibility$white, levels = c(1, 0), 
                                    labels = c("White", "Nonwhite"))
summary_credibility$treatment <- factor(summary_credibility$treatment, levels = c(0, 1), 
                                        labels = c("Control", "Treatment"))
p_cred_race <- 
  ggplot(data = summary_credibility, aes(x = white, y = cred, fill = treatment)) +
  geom_bar(stat = "identity", position = position_dodge(), color = "black", alpha = .5) +
  scale_x_discrete(breaks = c("White", "Nonwhite")) +
  scale_fill_manual(values = c("grey15", "grey85")) +
  geom_errorbar(width = .2, aes(ymin = conf.low, ymax = conf.high), 
                position = position_dodge(.9)) +
  xlab("") + 
  ylab("Perceived Credibility of Information\n(0 = very untrustworthy; 4 = very trustworthy)") +
  ggtitle("Perceived Credibility by Racial Identity") +
  theme_classic() +
  theme(text = element_text(family = "Times", size = 13),
        axis.text = element_text(color = "black", size = 13),
        legend.justification = c(1, 1), legend.position = c(.34, .99),
        legend.box.background = element_rect(color = "black"), 
        legend.key.size = unit(1, "line"),
        legend.direction = "vertical",
        legend.margin = margin(t = 0.2, r = 0.2, b = 0.2, l = 0.2, unit = "cm"),
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 13, face = "italic")) +
  coord_cartesian(ylim = c(0, 4))
credibility <- plot_grid(p_cred_pid, p_cred_race, labels = "AUTO", nrow = 1, label_fontfamily = "Times")
ggsave(file = "Figure S21.pdf", credibility, width = 9, height = 4.5)

### Figure S22 ----
ATE_pop_black <- lm_robust(pop_black ~ treatment + age + female + white + 
                             married + education + income + ideo + pid,
                           data = df)
group_mean_pop_black <- df %>% 
  group_by(treatment) %>% 
  summarize(xvalue = mean(pop_black, na.rm = T))
p1 <- 
  ggplot(data = df, aes(x = pop_black, group = treatment, fill = treatment)) +
  geom_density(alpha = .25) +
  geom_vline(data = group_mean_pop_black,
             aes(xintercept = xvalue, color = treatment), 
             linetype = "dashed", linewidth = 0.6) +
  scale_color_manual(labels = c("Control", "Treatment"), values = c("#00008B", "#EFC000FF")) +
  scale_fill_manual(labels = c("Control", "Treatment"), values = c("#00008B", "#EFC000FF")) +
  geom_vline(xintercept = 14.2, color = "red", linetype = "solid", linewidth = 1) +
  xlab("Estimated % of Blacks in the US Population") +
  ylab("Density") +
  ggtitle("Posterior Beliefs about Black Population in Full Sample") +
  theme_classic() +
  theme(text = element_text(family = "Times", size = 13),
        axis.text = element_text(color = "black", size = 13),
        legend.justification = c(1, 1), legend.position = c(.95, .99),
        legend.box.background = element_rect(color = "black"), 
        legend.key.size = unit(1, "line"),
        legend.direction = "vertical",
        legend.margin = margin(t = 0.2, r = 0.2, b = 0.2, l = 0.2, unit = "cm"),
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 13, face = "italic"))
p1 <- p1 + 
  annotate(geom = "text", x = 50, y = 0.032, 
           label = paste("beta == ", round(ATE_pop_black$coefficients[2], 2)),
           size = 4, family = "Times", parse = T) +
  annotate(geom = "text", x = 50, y = 0.030, 
           label = paste("(SE == ", round(ATE_pop_black$std.error[2], 2), ")"),
           size = 4, family = "Times", parse = T)
ggsave(file = "Figure S22.pdf", p1, width = 6, height = 4.5)

### Figure S20 ----
## ATE on affective polarization
group_mean_temp1 <- df %>% 
  group_by(group) %>% 
  summarize(xvalue = mean(aff_pol, na.rm = T))
ATE_aff_pol <- lm_robust(aff_pol ~ treatment + age + female + white + 
                           married + education + income + ideo + pid,
                         data = df)

## ATE on outparty affect
group_mean_temp2 <- df %>% 
  group_by(group) %>% 
  summarize(xvalue = mean(out_affect, na.rm = T))
ATE_out_affect <- lm_robust(out_affect ~ treatment + age + female + white + 
                              married + education + income + ideo + pid,
                            data = df)

## Visualize the distributions of affective polarization and outparty affect
p1 <- 
  ggplot(data = df, aes(x = aff_pol, group = group, fill = group)) +
  geom_density(alpha = .25) +
  geom_vline(data = group_mean_temp1,
             aes(xintercept = xvalue, color = group), 
             linetype = "dashed", size = 0.6) +
  scale_color_manual(labels = c("Control", "Treatment"), values = c("#00008B", "#EFC000FF")) +
  scale_fill_manual(labels = c("Control", "Treatment"), values = c("#00008B", "#EFC000FF")) +
  xlab("Affective Polarization\n(0 = lowest; 1 = highest)") +
  ylab("Density") +
  ggtitle("Affective Polarization") +
  theme_classic() +
  theme(text = element_text(family = "Times", size = 13),
        axis.text = element_text(color = "black", size = 13),
        legend.justification = c(1, 1), legend.position = c(.3, .99),
        legend.box.background = element_rect(color = "black"), 
        legend.key.size = unit(1, "line"),
        legend.direction = "vertical",
        legend.margin = margin(t = 0.2, r = 0.2, b = 0.2, l = 0.2, unit = "cm"),
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 13, face = "italic"))
p1 <- p1 + 
  annotate(geom = "text", x = 0.6, y = 1.75, 
           label = paste("beta == ", round(ATE_aff_pol$coefficients[2], 2)),
           size = 4, family = "Times", parse = T) +
  annotate(geom = "text", x = 0.6, y = 1.65, 
           label = paste("(SE == ", round(ATE_aff_pol$std.error[2], 2), ")"),
           size = 4, family = "Times", parse = T)

p2 <- 
  ggplot(data = df, aes(x = out_affect, group = group, fill = group)) +
  geom_density(alpha = .25) +
  geom_vline(data = group_mean_temp2,
             aes(xintercept = xvalue, color = group), 
             linetype = "dashed", size = 0.6) +
  scale_color_manual(labels = c("Control", "Treatment"), values = c("#00008B", "#EFC000FF")) +
  scale_fill_manual(labels = c("Control", "Treatment"), values = c("#00008B", "#EFC000FF")) +
  xlab("Outparty Affect\n(0 = lowest; 1 = highest)") +
  ylab("Density") +
  ggtitle("Outparty Affect") +
  theme_classic() +
  theme(text = element_text(family = "Times", size = 13),
        axis.text = element_text(color = "black", size = 13),
        legend.justification = c(1, 1), legend.position = c(.95, .99),
        legend.box.background = element_rect(color = "black"), 
        legend.key.size = unit(1, "line"),
        legend.direction = "vertical",
        legend.margin = margin(t = 0.2, r = 0.2, b = 0.2, l = 0.2, unit = "cm"),
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 13, face = "italic"))
p2 <- p2 + 
  annotate(geom = "text", x = 0.375, y = 2.3, 
           label = paste("beta == ", round(ATE_out_affect$coefficients[2], 2)),
           size = 4, family = "Times", parse = T) +
  annotate(geom = "text", x = 0.375, y = 2.15, 
           label = paste("(SE == ", round(ATE_out_affect$std.error[2], 2), ")"),
           size = 4, family = "Times", parse = T)
ATE_polarization <- plot_grid(p1, p2, labels = "AUTO", nrow = 1, label_fontfamily = "Times")
ggsave(file = "Figure S20.pdf", ATE_polarization, width = 9, height = 4.5)

### Figure S12 ----
## HTE on work ethic stereotype
out <- interflex(Y = "stereo", D = "treatment", X = "overest",
                 Z = c("age", "female", "white", "married", "education", "income", "ideo", "pid"), 
                 data = df, na.rm = TRUE,
                 estimator = "kernel", vcov.type = "robust")
p1_HTE <- 
  plot(out, theme.bw = TRUE, cex.axis = 0.8, cex.lab = 0.8, bin.labs = TRUE,
       xlab = "Level of Overestimation", 
       ylab = "Marginal Effect on Work Ethic Stereotypes",
       ylim = c(-1, 1)) +
  ggtitle("HTE on Work Ethic Stereotypes by Level of Overestimation") + 
  theme(plot.title = element_text(hjust = 0.85, size = 13, face = "italic", family = "Times"))

## HTE on racial resentment
out <- interflex(Y = "resent", D = "treatment", X = "overest",
                 Z = c("age", "female", "white", "married", "education", "income", "ideo", "pid"), 
                 data = df, na.rm = TRUE,
                 estimator = "kernel", vcov.type = "robust")
p2_HTE <- 
  plot(out, theme.bw = TRUE, cex.axis = 0.8, cex.lab = 0.8, bin.labs = TRUE,
       xlab = "Level of Overestimation", 
       ylab = "Marginal Effect on Racial Resentment",
       ylim = c(-1, 1)) +
  ggtitle("HTE on Racial Resentment by Level of Overestimation") + 
  theme(plot.title = element_text(hjust = 0.85, size = 13, face = "italic", family = "Times"))
HTE_overest <- plot_grid(p1_HTE, p2_HTE, labels = "AUTO", nrow = 1, label_fontfamily = "Times")
ggsave(file = "Figure S12.pdf", HTE_overest, width = 10, height = 4.5)

### Table S8 ----
## CACE on work ethic stereotype
cace_1 <- ivreg(stereo ~ mp_check | treatment, data = df)
cace_2 <- ivreg(stereo ~ mp_check + age + female + white + married +
                  education + income + ideo + pid | 
                  treatment + age + female + white + married +
                  education + income + ideo + pid, data = df)

## CACE on racial resentment
cace_3 <- ivreg(resent ~ mp_check | treatment, data = df)
cace_4 <- ivreg(resent ~ mp_check + age + female + white + married +
                  education + income + ideo + pid | 
                  treatment + age + female + white + married +
                  education + income + ideo + pid, data = df)

## Summarize the results
m_list <- list(IV1 = cace_1, IV2 = cace_2, IV3 = cace_3, IV4 = cace_4)
msummary(m_list, output = "markdown")
